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Abstract 

Collective excitations in simple metal systems can be described successfully in 
terms of a local one-body excitation operator Q, due to the long range nature of 
the coulomb interaction. For the plasmon modes of a simple-metal slab, momen- 
tum expansions of Q are calculated using a variational procedure, equivalent to a 
restricted RPA calculation. The dispersion relation and the density fluctuation for 
each mode are found in the sudden approximation using the proper Q operator and 
the RPA sum-rule formalism. The contribution of the exchange and correlation 
energy is estimated using a local density functional. The positive background is 
described within a jellium model while the ground-state electronic density is ap- 
proximated by a double step profile. The density fluctuation of the plasmon modes 
above the plasma frequency form standing waves across the slab. The spectra below 
the plasma frequency is qualitatively different to that of local optics calculations, 
due to the appearance of two multipole plasmon modes that shift down the origin 
of the ui+ plasmon. The dependence of the results on the width of the slab, the 
density of the simple-metal and the surface diffuseness is discussed. Throughout, 
we compare our results with previous RPA and TDLDA calculations. 
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1 Introduction. 



There are several ways in which the dispersion relation of the plasmon modes in systems 
with translational invariance can be measured. The first experiments on simple metal 
slabs used a dielectric substrate to support the metal film and measured the light 
transmittance pj or the photoelectric yield ||. More recently, molecular-beam epitaxy 
has been applied to grow thin metallic layers between two insulators Q or Al x Ga\- x As 
heterostructures simulating a positive jellium background Also electron energy 
loss spectroscopy (EELS) has been used to study adsorbed simple-metal layers on an 
aluminium substrate ||. 

The spectrum of a simple-metal slab has two distinct parts: the plasmon modes 
that lie above ui p , corresponding to standing waves across the slab, and the surface 
plasmon modes, with lower frequencies. Some theoretical calculations have been done 
using the random-phase approximation (RPA) || [/], ||, |9|, |10| and the time-dependent 
local density approximation (TDLDA) [|ll], [l^, [l3|]. The lo- and the plasmons are 
studied thoroughly in all this references, but only two of them discuss the existence of 



an additional multipole plasmon mode. Liebsch [12], in a calculation of the plasmon 
modes of a simple-metal overlayer on an aluminium substrate, finds a multipole plasmon 
that hybridizes with the u>+ plasmon. Schaich and Dobson |l3| also obtain a peak in 
the response function at frequencies between u> p and u> p / v2 but they do not conclude 
whether the cause is a collective mode or is due to intersubband excitations. 

This multipole plasmon mode has been detected in simple-metal clean surfaces 



using EELS [14]. Its dispersion relation has a positive slope when the momentum k 
parallel to the surface is zero while the frequency at this point is 0.8ui p . This value 
agrees with theoretical calculations 14, [15|]. Multipole plasmon modes are associated 
with electronic-density fluctuations that are peaked at the surface region and have 
decreasing oscillating amplitude towards the interior of the metal. It can be shown [16] 
that the integral of the electronic density perpendicular to the surface is zero. Unlike 
ordinary surface plasmons, they are optically active because they carry momentum on 
the normal direction. All this properties were found for simple-metal surfaces but one 
expects that at least some of them should be valid for slabs. 

The aim of this work is to study the existence of multipole plasmon modes on simple- 
metal slabs and its coupling to the other surface modes. We will use the RPA sum-rule 
formalism (SR) that previously has been able to produce useful results for different 
geometries U, HI, |19|, ||. In particular, a good agreement has been obtained between 
experimental data and the dispersion relation of surface plasmons [21] and multipole 
plasmons [15] on a plane simple- metal surface. It has some advantages over other 
methods. The different contributions to the energy of each mode can be analyzed on 
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its own and semiclassical approximations for the kinetic energy are easily implemented. 
It can incorporate electronic exchange and correlation interactions using a local density 
approximation. Within the SR formalism it is not difficult to employ realistic electronic 
profiles. In fact the results have RPA precision as long as a self-consistent ground-state 



electronic profile is used [22], although only an average of the response function can be 
obtained. 

Usually, when working within the SR method, it is assumed that for every collective 
mode exists a local one-body observable that excites that mode without mixing it with 
others. In order to obtain useful results an appropriate excitation operator Q must 
be found and the system must have well-defined collective modes. Otherwise it would 
not be possible to get any valuable information from the averaged response function. 
There is not a unique way to find such an operator. Sometimes the properties of a given 
operator Q determine the kind of mode that it excites. For example, an operator that 
fulfills Laplace's equation can only excite surface modes pi]] . A restricted variation of 
Q on a parameter [15 can be useful if a sufficiently general form of Q is employed. In 
RPA calculations an excitation operator Q can be found for each mode but although it 
will be an observable, in general it is not a local operator. Nevertheless, the fact that 
the interaction responsible for the existence of well-defined collective excitations, the 
coulomb interaction, has a long range, allows to approximate Q as a local operator. This 
approximation is called local- RPA 23] and has been applied successfully to nuclei [24] 
and metal clusters We use local- RPA to find a proper Q operator as a momentum 
expansion. 

Instead of a self-consistent electronic profile, a double step profile is used to de- 
scribe the ground-state. This leads to a simplification of the calculations which can 
be performed analytically without a great loss of precision because it has been shown 
[21] that there is little quantitative difference between the results from any of the two 
models. 



2 Application of local RPA to a metal slab 

In this section we first give a short review of the method as it was stated in [ p3| ] (see 
this reference for more details). Then we expatiate on the particularities that we have 
had to devise to adapt it to our problem, the most important being the choice of the 
form of the functions that give rise to the sought-after excitation operators. 

Both local and full RPA demand as a starting point a description of the vacuum 
of the system at issue. Within the quasi-bosonic approximation this description is fur- 
nished by working out the Hartree-Fock approximation to the vacuum state \(f>o). This 
approximation is acceptable as long as one assumes that correlations are small. Once 
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the Hartree-Fock problem is solved, the essence of the method is to build the excitation 
operators as linear combinations of one-particle-one-hole operators, and determine the 
coefficients of the expansion by solving the equation of motion 

(00 1 C m ,[H,Cl] |0 O ) = S nm hu m , (1) 

where C\ is an excitation operator. By this name we mean that when it acts on the 
ground state, an elementary excitation of the system with a well defined energy is 
created. An alternative description of RPA can be made based upon the operator 

Q,-^ (2) 

H,[H,Q n ] ph ] ph = (Hu n ) 2 Q n , (3) 

The ph symbol means that the operator is projected onto the linear space spanned by 
one-particle-one-hole states. Local RPA assumes that the operator Q n can be expanded 
in terms of a function /„ as follows, 

N 

Qn = £/n(n), (4) 

1=1 

where N is the number of electrons and = (xj , yi , zi) is the position operator of the 
i-th particle. In order to determine f n we shall consider that it is a linear combination 
of the elements of a basis of functions that we shall choose taking advantage of the 
particular problem we are facing at. We will focus on the problem of a gas of electrons 
in a simple metal slab of very large area A and thickness d. If the gas is made of N 
electrons, then the operator of total momentum q to t of the gas will be; 

N 

qtot = X) & (5) 
i=i 

We will assume the slab to be parallel to the XY plane and consequently write any 
vector as q = (k,p), k being the 2 dimensional vector that lives in the x — y plane and 
p is the projection along the z axis. Translation invariance in the x and y axis leads to 
the conservation law; 

H, k tot ] = (6) 

So, the eigenvalues k of operator k are good quantum numbers for our system. As k is 
a true quantum number of the problem, we may advance that each Q n operator has a 
well defined k, in other words, it has the form 

N 

Q n = J2e^M^)- (7) 
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This observation suggests that we relabel Q n as Qk,n hi the sequel. 

As for the /„ functions, which we haste to relabel in accordance with the previous 
discussion, we may write 

fk, n (R, z) = exp (ikR)(j) n (z). (8) 

The electronic profile vanishes for distances far enough from the boundaries of the slab. 
We can think the gas of electrons to live in a linear box of size L. We must take L to 
be larger than d as there are electrons wandering outside the slab, but provided L is 
large enough, the physical results are L-independent. Nevertheless the value of L must 
not be larger than necessary because this would increase the number of terms needed 
in the expansion. 

Now we have a physically motivated choice for the basis which generates <p n ; we 
take cj) n to be expanded in terms of a Fourier basis 

<t> n tz) = Y,a l nj piz , (9) 
i 

with pi = ^zl. Equivalently, we may expand the Qk,n operators 

Qk,n = ^2a l n Qk, Pl , (10) 
i 

with 

N 

Qk,Pi = ^2 ex P ( ik Rj) ex P (ipizj)- ( n ) 

By construction, the Q operator will have a period L. This is done just for mathematical 
convenience. We have truncated the previous Fourier series ((I) or ©) making sure 
that the terms we do not consider are really negligible. This can be done checking the 
convergence of the results as the number of terms is increased. It is clear that L cannot 
be too small as we would cut the tails of the electronic density, but it should not be 
larger than necessary because this would increase the number of terms needed in the 
preceding expansion, a? are coefficients that can be obtained by minimizing the energy 
functional [23]. Written in terms of the Fourier coefficients eq. |3| reads 



K a p ~ (hLo n ) 2 B a/3 \ < = 0, (12) 
JC a p and B a/ 3 being respectively; 

B^ = (^ \[Q a ,[H,Qp}}\%}, (13) 

tC a p = <$o| [[Qa, H] , [H, [H, Qp]]] |$o>- (14) 
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H is the Hamiltonian of the system and |$o) the (Hartree-Fock) ground state. Needless 
to say that the sum-rule approach 



m 3 



m 1 = h$ \[Q,[H,Q])\$t 



^(%\[[Q,H],[H, [H,Q]]] |$ ), 



(15) 

(16) 
(17) 



is a particular case of eq. 12 because if we could find a series expansion such that every 



term would couple to a different eigenmode, eq. 12 would be separated into a set of 



equations like eq. 17, one for each term of the series. It can be shown [27] that E% is 
an upper bound to the energy of the mode excited by Q. 



Solving equation |12| for a fixed value of k does not only deliver the optimal operator 
Q but also the dispersion relations of the collective excitations w n (k)- We calculate the 
dispersion relation in an equivalent, but more transparent way. As the Hamiltonian 

-,2 



H 



dr 



7r 



2m, 
= Hk + H c . 



, . n(r) ( 9 f , , n(r') 
- T ( r ) + -pr 2 ( e / dr ~ ,,, + 2 Vi (r 



(18) 



contains a kinetic and a coulomb term, we can evaluate separately the contribution of 
each part of the Hamiltonian to the 771,3 sum-rule. When the Q operator used is that 
obtained solving eq. Il2|, the eigenvalues of this equation are, by construction, the £3 



energy obtained from eq. 17 taking into account all the contributions to 7773. If E3 is 
calculated using only a contribution to 7773, the quantity obtained gives an idea of the 
relative importance of that term in the total energy. The inclusion of an exchange and 



correlation term when solving eq. 12 is done in a perturbative way, as explained later 
on. 

Another important magnitude in the sum-rule formalism is the sudden approxima- 
tion to the density fluctuation given by 



fiV nVQ 



(19) 



It is the first-order contribution to the density fluctuation of the |7y) state described in 
the appendix and usually gives a reliable description of the eigenstates of the system. 



At this point we just need to solve equation 12 and we will do so by expressing JC a p 
and B a j3 in terms of the local density n(r) and the kinetic-energy density r(r) (due to 
the symmetry discussed so far the only dependence of this quantities is on the z axis) 



t r 



(20) 
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To compute the kinetic term the most direct approach consists in taking the kinetic 
energy of a free gas of electrons; 

T = ara(r) 5 / 3 . (21) 

The coulomb part contains the coulomb interaction between electrons (of charge e) 
and the interaction with the jellium through a potential v(r). An estimation of the 
exchange and correlation energy is made using the following Slater- and Wigner-type 
expressions 

e ex = —an{r) 1 ^ - b \/ 22 
4 cn(r) i /' 3 + a 



with a = y/3/n,b = A4,c = 7.8 and d = a/2. Although this term is not included in the 
Hamiltonian when computing Q, a contribution to the to 3 sum-rule is calculated just 
as if it were another part of eq. O and added to the other contributions. So the final 



values of the dispersion relation are slightly different to the ones obtained in eq. 12 
In this way the exchange and correlation interactions are included in the results. The 
main drawback is that the self-consistency of the calculation is lost to a certain extent 
but this is acceptable because exchange and correlation do not play a crucial role in 
the existence of the collective excitations in an electron gas. 

The details necessary to solve equation ^ and improvements to the kinetic energy 
functional are postponed to the appendix as they become rather technical and do not 
contain any new physics. 

3 Results and discussion 

We have applied the method described in the preceding section to simple metal slabs. 
A double step electronic profile is used as a description of the ground state 



n[z) = — 
V ' 2 



n, d + 5. „, d — 5. d + 5. n . d — 5, 
9(z + — — ) + 9{z + — — ) - 9(z — ) - 9{z — 



(23) 



The length of the step 5 has been adjusted in all the cases by a minimal squares pro- 
cedure to an improved Thomas-Fermi electronic profile |2(| calculated for each system. 
In this section atomic units are used throughout. 

Figure 1 (a) represents the dispersion relation of the plasmon modes below lo p for 
a typical system. The lower mode is the oj- mode. It begins at w = for k = 
and tends to the surface plasmon frequency Up/y/2 as k tends to infinite. This 
behaviour is completely in agreement with Local Optics, not like the next mode, that 
we identify as the iv + mode of Local Optics. In our calculation it has a frequency of 
approximately to = 0.8w p for k = when according to Local Optics, and most RPA 
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and TDLDA calculations, it should tend to uj = uo p for a small k. For large values of 
k it becomes degenerate with the ui- mode. There is another mode whose frequency 
tends to 0.8co p for a low k. We call it the even multipolar mode and it has a qualitative 
behaviour very similar to that of multipolar modes in a surface. Its dispersion relation 
has a positive slope when k = and tends to the same value as in the semi-infinite 
case. In our calculation it becomes degenerate for large k with a mode that we call 
the odd multipolar mode. This last mode has its origin at the plasma frequency uj p . 
Although this scheme is not the usual one in semiclassical or fully quantum-mechanical 
calculations of collective excitations in metal slabs, a nearly identical behaviour has 
been reported for plasmons in a simple-metal overlayer [12|. Figure 1 (b) displays 
the dispersion relation of the first four bulk resonances. These dispersion relations 
are qualitatively similar to those obtained by discretization of the z component of the 
momentum in a bulk plasmon dispersion relation but the slopes are flatter than those 
obtained with other methods, even classical ||. Also the position of the origins does 
not seem to be quadratic as a reasoning of this type would indicate. 

Further insight can be gained looking at the sudden approximation of the density 
fluctuation for each mode. Every one of the four discontinuities in the ground-state 
electronic profile produces an avoidable discontinuity in the density fluctuation. In 
those points the first derivative of the density fluctuation is also discontinue. A more 
realistic profile should be used to obtain quantitatively correct results but as the double- 
step model includes the most important physical aspects, such as diffuseness, qualitative 
information like symmetry and number of nodes is reliable. The density fluctuation of 
the uj- and the u;+ modes is shown in Fig. 2 (a). The symmetry and the number of 
nodes of this modes are the same as those of the even multipole and the odd multipole 
respectively, in Fig. 2 (b). All these modes are clearly surface modes and the main 
difference between them is that in the multipole modes the density fluctuation is not 
constant in the middle of the slab. The decreasing oscillating amplitude that multipole 
plasmons show in a semi-infinite medium is not found in our calculation for thin slabs. 
The density fluctuation of the bulk resonances is represented in Fig. 2 (c)-2 (d). The 
first resonance has two nodes and presents the typical profile of a standing wave across 
the slab. In fact, we can label a resonance with its number of nodes because any 
resonance has one more node than the previous one. 

The key part of the calculation is to choose the number M of terms of the Q ex- 
pansion. This is done, for a given momentum k parallel to the surface, examining the 
energy of the modes in front of M. The value of M necessary to attain the convergence 
depends on a great number of physical factors. As the main assumption of Local RPA 
is the existence of well defined collective excitations, the convergence will be difficult or 
even impossible in all those situations where a severe damping of the plasmon modes 
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is expected. Because the relative importance of the coulomb energy increases with 
the electronic density, in high-density metals like aluminium the convergence will need 
lower values of M than low-density metals like cesium. Also when the momentum k 
approaches the Landau cutoff value, given approximately by q c = u p /vf, the conver- 
gence must be harder. The surface diffuseness also plays its role because an increasing 
diffuseness tends to lower the coulomb energy. All these general trends are present in 
our calculations. The size of the system is also an important factor. This is mainly 
because the expansion series of a bigger system needs more terms to be accurate. This 
puts un upper bound on the size of the slabs that can be treated within this method. 

Taking advantage of the transparency of the method, the contributions of the kinetic 
and the coulomb parts of the hamiltonian to the dispersion relation are represented in 
Fig. 3. The ratio between the coulomb and the kinetic energy is different for each mode. 
The mode has a negligible kinetic energy in front of the coulomb contribution while 
they are comparable in the cu+ mode. In fact, looking at Fig. 3 one can expect the w+ 
mode to have the worse convergence of all, as it really happens in our calculation. This 
agrees with TDLDA results that show how the peaks in the excitation spectra are 
noticeably wider for the w+ mode than for the w_ mode. 

The dependence of the dispersion relation on the width d of the slab can be seen 
in Fig. 4 (a)-4 (b). The lo- and the to + modes tend to degenerate when k » d^ 1 and 
the same is true for the even and the odd multipole modes. This is due to the fact that 
the symmetry of the mode is not a real difference when the density fluctuation of two 
modes is concentrated on the two surfaces of the slab and these are enough apart in 
comparison with the wavelength of the mode. Obviously, for larger values of d plasmons 
will tend to degenerate for smaller values of k. The lo- and the a>+ modes degenerate at 
an energy close to u p /y/2. In the semi-infinite medium limit their dispersion relations 
will join at the origin and they will form the well known surface plasmon. The even 



and the odd multipole modes will also join in a unique multipole plasmon mode [15| for 
extremely large slabs. In Fig. 4 (b) can be seen that the dispersion relation of the odd 
plasmon mode is much more affected by the change of the width of the slab than that 
of the even multipole mode, so the dispersion relation of the surface multipole mode 
will be closer to that of this last mode. 

In Fig. 5 (a) the dispersion relations of two slabs with a different electronic density 
are plotted. The sodium slab has a greater electronic density and the energy of its 
modes is higher than those of the potassium slab. Surface modes degenerate for lower 
values of k in slabs of this last metal. The increase in the energy of the modes in a 
high-density metal is due to both kinetic and coulomb energies although in relative 
terms the coulomb term is more responsible of this feature. Moreover, the density 
also influences the diffuseness of the electronic profile. As discussed earlier, these two 
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factors have an important effect on the convergence of the expansion series. A plot of 
the plasmon modes for an aluminium slab is displayed in Fig. 5 (b). Besides the rise 
of the energies of the modes we can notice a greater gap between the origin of the odd 
multipole mode and the first bulk resonance. This is not due directly to the higher 
electronic density but to the steeper slope of the electronic profile as will be seen soon. 
In order to distinguish the effects due to the change of electronic density and those 
caused by a different diffuseness of the profile we have studied three systems where we 
have changed only the length of the intermediate step. The results are represented in 
Fig. 6. When the length of the step is reduced the origins of the u>+ and of the even 
multipole tend to rise, they are no longer at 0.8w p , and a gap appears between them. 
It is clear that the origin 0.8w p depends strongly on the diffuseness of the surface and 
can be used as a test on how good is an electronic profile. The odd multipole and 
the first bulk resonance also have higher energies for k = and the gap between their 
origins also increases. The slopes of the dispersion relations of the even and the odd 
multipole modes suffer a significant change and are closer than before. Finally, for a 
system with no diffuseness, the origin of the w+ mode has risen up to oj v and the two 
dispersion relations of the multipole modes are in the same place as the first and second 
bulk resonances in the previous systems. The result of the reduction of the diffuseness 
is, then, the gradual conversion of the even and the odd multipole modes into the first 
two bulk resonances. 



4 Summary and conclusions 

We have succeeded in developing a semi-classical calculation that includes as solutions 
the even and the odd multipole modes, that are not usually considered in related 
calculations. In the semi-infinite medium they are brought together into the by now well 
studied multipole plasmon mode. These multipole modes modify the classical origin 
of the uj+ mode and are closely dependent on the diffuseness of the electronic profile. 
Their origin is shifted upwards when the diffuseness is reduced and they disappear 
when it is set to zero. 

The u>- plasmon is far more observable, for intermediate values of k, than the u>+ 
plasmon. This can be stated either examining the convergence of the results on the 
number of terms of the expansion series or taking into account the relative contribution 
of the coulomb and the kinetic energies. The bulk resonances closer to uj p are also 
realistically described. The Thomas-Fermi density functional gives a too flat slope of 
the kinetic energy. The exchange and correlation contribution lowers the energy for 
large values of k. 

Local RPA is a useful and economic method for finding the excitation operator 
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Q whenever well-defined collective excitations are expected. The size of the system 
determines, among other factors, the number of terms needed in the expansion series. 
If a strong damping acts on the plasmon modes due to high diffuseness of the electronic 
profile, low electronic density or whatever other cause, the convergence can be difficult 
to reach or even non-existent at all. 
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A Appendix 

In this appendix we will present the formulae which were used in the solution of equation 
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The first object to compute, that is £>, is a straightforward calculation; 

B(k,p),(k>, P >) = M k2 + PP') I dz exp (-i(p - p')z)n(z) (24) 



To compute K, a p we start by considering; 

M al3 = M a/3 k + M a(i ex + M af f 

= ($ | exp V ([Q a ,H]) exp (7/ [Q^, #]) (H k + F ei . + ff c ) x 

exp (-7/ [Qp,H]) exp (-77 [Q a ,H])\$ ) 
= (rjri'lHk + H^ + H^rjrj') (25) 

The previous object is useful as it is connected to the one we are interested in through 
formulae; 

JC a p = ( JjnMafs ) (26) 



I drjdr]' " 



v =0, v '=0 



At this point we have reduced the problem of calculating the expressions in for- 
mula 14 to the evaluation of expressions such as {r]rj'\F(n)\r]rj') which are still quite 
intractable. Some approximations are needed to go any further. We will approximate 
terms like this one as 

(777/ 1 1 7777') = F( (7777' 1 71)7777')) + ... (27) 

and keep just the first term. 

The rest of the computation is straightforward but lengthy. We arrive at 

= ^App' / dzn(zf exp(-i(p - p')z) 

^/c 3 — pp'k + i(p + p')k 2 sgn(z — z')J 

—AirApp J dz (n(z) — rij(z')) n(z) exp(—i(p — p')z) 

-2i,A (pp> + #) p> J dzdz* (n(z) - nj ( Z ')) n(z) exp(-^ - p>)z)s g n(z 
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The exchange term is easily computed with the result; 

-(k 2 + p' 2 )p{p - p') J dzexp(-i(p-p')z) ( n ( z ) ^^y ~ e x{n)j 

(29) 

The kinetic part does not show up any further complication as compared to the ex- 
change term since where one reads e ex one substitutes r(n) and the rest is exactly the 
same. 
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Figure captions 

Figure 1 (a): Dispersion relation of the plasmon modes below the plasma frequency 
in a sodium slab [d = 25ao where ao is the Bohr radius). From bottom to top can be 
seen the u—, u;+,even multipole and odd multipole modes. Thin lines represent, from 
bottom to top, w p /\/2, 0.8oj p and u p . 

Figure 1 (b): Dispersion relation of the four lowest bulk plasmon resonances in 
the same case as Fig. 1 (a). The thin line represents the plasma frequency. 

Figure 2 (a): Sudden approximation to the density fluctuation calculated as ex- 
plained in the text for a sodium slab (d = 40ao). The dashed line represents the w_ 
mode while the solid line is the u> + mode. 

Figure 2 (b): The dashed line represents the density fluctuation of the even mul- 
tipole mode in the same case as Fig. 2 (a). The solid line is the odd multipole mode. 

Figure 2 (c): Density fluctuation of the two lowest bulk plasmon resonances for 
the same system as Fig. 2 (a). The dashed line represents the lowest resonance and 
the solid line represents the next resonance. 

Figure 2 (d): The dashed line is the third bulk resonance with lowest energy and 
the solid line is the fourth. Calculation for the same system as previous cases. 

Figure 3: The energy contributions of the same mode are plotted in the same 
type of line. The upper one is the coulomb contribution and the lower one is the 
kinetic contribution. The solid line corresponds to the w_ mode, the long-dashed to 
the cu+ mode, the short-dashed to the even multipole mode, the dotted line to the odd 
multipole mode and the dot-dashed line corresponds to the lowest bulk resonance. The 
calculations are made for the same case as Fig. 1 (a). 

Figure 4 (a): Dispersion relation of the uj- and the u;+ modes for a sodium slabs 
of different widths. The solid line represents a width 18ao, the long-dashed line stands 
for 25ao and the short-dashed line is the result for 40ao- 

Figure 4 (b): Same as Fig. 4 (a) but with the even and the odd multipole plasmon 
modes. 

Figure 5 (a): Dispersion relation of the plasmon modes of two slabs of different 
metals and the same width (d = 25ao). The solid line corresponds to a potassium slab. 
The dashed line is used for the sodium slab. 

Figure 5 (b): Dispersion relation of the plasmon modes of an aluminium slab. 
The width is 25ao- 

Figure 6: Dispersion relation of the four modes with lowest energy. The solid line 
represents the modes with an intermediate step set to S = 2.11ao- The long-dashed 
lines are calculated with 5 = Ioq while while short-dashed lines stand for a system with 
5 = 0. 
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